SoSe2021

Folienübersicht

Lineare Modelle in R

Anwendungsbeispiel iris

lm()

Wir modelliere die Kronblattlänge als Funktion der Kelchblattlänge:

mod <- lm(formula = Petal.Length ~ Sepal.Length, data = iris)
mod
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)

Coefficients:
 (Intercept)  Sepal.Length  
      -7.101         1.858  
coef(mod)
 (Intercept) Sepal.Length 
   -7.101443     1.858433 

Wie lautet die vollständige lineare Regressionsgleichung?

Populations- vs. Stichprobengleichung

1

\[y = a + b*x\]

  • Die häufigste Art, die Gleichung zu schreiben ist die Stichprobengleichung für die erklärenden y-Werte!



Zur Erinnerung

  • Parameter = eine Messgröße, die eine POPULATION beschreibt (z.B. \(\mu, \sigma^{2}\))
  • Statistik = eine Messgröße, die eine STICHPROBE (‘sample’) beschreibt (z.B. \(\bar{x}, s^2\))

Populations- vs. Stichprobengleichung

2

Populationsgleichung:

  • für beobachtete y-Werte: \(y_{i} = \alpha + \beta*x_{i} + \epsilon_{i}\)
  • für vorhergesagte y-Werte: \(\mu(y_{i}) = \alpha + \beta*x_{i}\)

Stichprobengleichung:

  • für beobachtete y-Werte: \(y_{i} = a + b*x_{i} + \epsilon_{i}\)
  • für vorhergesagte y-Werte: \(\hat{y}_i = a + b*x_{i}\)

Die Stichprobenkennwerte \(a\) und \(b\) repräsentieren die Punktschätzer der Regressionsparameter \(\alpha\) und \(\beta\)!

Anwendungsbeispiel iris

coef(mod)
 (Intercept) Sepal.Length 
   -7.101443     1.858433 

Für die beobachteten Stichprobenwerte gilt:

\[Petal.Length_i = -7.1 + 1.85*Sepal.Length_i + \epsilon_{i},~~\text{wobei}~~\epsilon_i ~N(\mu, \sigma^2)\]

R rechnet noch viel mehr aus: die summary() Funktion

summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.47747 -0.59072 -0.00668  0.60484  2.49512 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -7.10144    0.50666  -14.02   <2e-16 ***
Sepal.Length  1.85843    0.08586   21.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared:   0.76, Adjusted R-squared:  0.7583 
F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

Wo ist a? Wo ist b? Und was bedeuten die anderen Zahlen…?

Die Genauigkeit der Schätzung

Von der Stichprobe auf die Grundgesamtheit verallgemeinern

  • Wie wären die geschätzten Koeffizienten, wenn wir das Experiment wiederholen würden?
  • Wo liegen die geschätzten Koeffizienten für die gesamte Population?
  • Hier braucht es die Berechnung der Standardfehler und Hypothesentests
  • Zur Erinnerung: Die Genauigkeit der Punktschätzung sollte immer durch den Standardfehler des Schätzers angegeben werden → je kleiner der Standardschätzfehler, desto größer ist die Effizienz des Schätzers.

Reproduzierbarkeit geschätzter Koeffizienten

Hohe Varianz

Stelle Dir vor, jeder Datenpunkt repräsentiert einen individuellen Fisch und alle Datenpunkte zusammen repräsentieren eine Population.

Wenn Du jetzt mehrere Zufallsproben von 10 Fischen nimmst, wie ähnlich wären sich dann die geschätzten Koeffizienten?

Reproduzierbarkeit bei starker Streuung

N = 10

Rate: Welche sind sich viel ähnlicher und welche sind näher an den wahren Populationskoeffizienten dran?

Reproduzierbarkeit bei starker Streuung

N = 10

‘Wahre’ Parameter: \(\alpha\) = 11.27 und \(\beta\) = 0.41

Reproduzierbarkeit bei starker Streuung

N = 100

Lass uns jetzt die Probengröße auf 100 Fische steigern!

Reproduzierbarkeit bei starker Streuung

N = 100

‘Wahre’ Parameter: \(\alpha\) = 11.27 und \(\beta\) = 0.41

Reproduzierbarkeit bei schwacher Streuung

N = 10

Wie ähnlich sind sich die geschätzten Koeffizienten bei schwacher Streuung und einem Stichprobenumfang von 10 Fischen?

Reproduzierbarkeit bei schwacher Streuung

N = 10

‘Wahre’ Parameter: \(\alpha\) = 10.04 und \(\beta\) = 0.5

Reproduzierbarkeit bei schwacher Streuung

N = 100

Lass uns jetzt die Probengröße auf 100 Fische steigern!

Reproduzierbarkeit bei schwacher Streuung

N = 100

‘Wahre’ Parameter: \(\alpha\) = 10.04 und \(\beta\) = 0.5

Maß der Genauigkeit der geschätzten Koeffizienten: \(s_\bar{a}\) und \(s_\bar{b}\)

Die Unbestimmtheit der geschätzten Steigung und Achsenabschnitt

  • steigt mit zunehmender Varianz und abnehmender Probengröße
  • ist größer wenn die Spanne der Werte klein ist (gemessen als \(SS_{X}\))

Standardfehler des Achsenabschnitts und der Steigung:

\[s_\bar{a}=\sqrt{\frac{MS_{Residual}\sum x_i^2}{n*SS_x}}=\sqrt{\frac{\frac{SS_{Residual}}{df}\sum x_i^2}{n*SS_x}}=\sqrt{\frac{\frac{\sum(y_i-\hat{y_i})}{n-2}\sum x_i^2}{n*\sum(x_i-\bar{x})^2)}}\]


\[s_\bar{b}=\sqrt{\frac{MS_{Residual}}{SS_x}}=\sqrt{\frac{\frac{SS_{Residual}}{df}}{SS_x}}=\sqrt{\frac{\frac{\sum(y_i-\hat{y_i})}{n-2}}{\sum(x_i-\bar{x})^2}}\]

Zur Erinnerung: Varianzzerlegung in linearen Modellen

Die Gesamtstreuung (\(SS_Y\) oder \(SS_{Gesamt}\)) kann in 2 Komponenten zerlegt werden:

erklärbare Streuung = \(SS_{Regression}\) + nicht erklärbare Streuung = \(SS_{Residuen}\)

Was sind die Standardfehler im Kronblattmodell?

summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.47747 -0.59072 -0.00668  0.60484  2.49512 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -7.10144    0.50666  -14.02   <2e-16 ***
Sepal.Length  1.85843    0.08586   21.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared:   0.76, Adjusted R-squared:  0.7583 
F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

→ Punktschätzung für den Achsenabschnitt: -7.10144 ± 0.50666

→ Punktschätzung für die Steigung: 1.85843 ± 0.08586$

Konfidenzintervalle der ‘wahren’ Koeffizienten

confint()

  • Konfidenzintervalle kann man nicht nur für Mittelwerte oder Varianzen berechnen, auch für Regressionskoeffizienten.
  • Berechnung basiert auf der t-Verteilung:

Manuelle Berechnung

a <- -7.10144
b <- 1.85843
ci_a <- 0.50666 * qt(p = 0.975, df=148)
ci_b <- 0.08586 * qt(p = 0.975, df=148)
mat <- matrix(
  c(a-ci_a, a+ci_a, b-ci_b, b+ci_b), 
  nrow = 2, byrow = TRUE)
rownames(mat) <- c("(Intercept)", "Sepal.Length")
colnames(mat) <- c("2.5%", "97.5%")
mat
                  2.5%     97.5%
(Intercept)  -8.102662 -6.100218
Sepal.Length  1.688760  2.028100

Mit build-in Funktion

confint(mod, level = 0.95)
                 2.5 %    97.5 %
(Intercept)  -8.102670 -6.100217
Sepal.Length  1.688772  2.028094
  • Mit 95%iger Wahrscheinlichkeit liegt \(\alpha\) zwischen -8.1 und -6.1 und \(\beta\) zwischen 1.7 und 2.0
  • Beide KI schließen nicht die Null mit ein → wir können davon ausgehen, \(\alpha \neq 0\) und \(\beta \neq 0\).

Signifikanztest der Koeffizienten

  • Die Schätzungen der Koeffizienten basierend auf der Stichprobe weichen immer von Null ab.
  • Ob dies Zufall ist oder nicht kann mit einem t-Test jeweils geprüft werden:
    • \(H_0:~\alpha = 0\), d.h. die Regressionsgerade geht durch den Ursprung.
    • \(H_0:~\beta = 0\), d.h. der Mittelwert ist für die Schätzung ebenso geeignet wie die Regressionsgerade.
  • Die Teststatistik t ist einfach nur der Quotient aus dem Schätzwert und seinem Standardfehler und folgt einer t-Verteilung:
    • \(t_a = \frac{a}{s_\bar{a}}~~\text{und}~~t_b=\frac{b}{s_\bar{b}}\)

Teststatistiken

(t_a <- -7.10144 / 0.50666)
[1] -14.01618
(t_b <- 1.85843 / 0.08586)
[1] 21.64489

Wahrscheinlichkeiten

pt(t_a, df = 148)
[1] 3.065741e-29
pt(t_b, df = 148, 
  lower.tail = FALSE)
[1] 5.224047e-48

Testergebnisse im summary() output

summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.47747 -0.59072 -0.00668  0.60484  2.49512 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -7.10144    0.50666  -14.02   <2e-16 ***
Sepal.Length  1.85843    0.08586   21.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared:   0.76, Adjusted R-squared:  0.7583 
F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

Interpretation der p-Werte

Wenn ein p-Wert sehr niedrig ist ( < 0.05), bedeutet das, dass

  1. der wahre \(\alpha\) und \(\beta\) Wert nicht 0 ist.
  2. Du eine ungewöhnliche Stichprobe hast.
  3. Deine Annahmen für das Modell falsch sind → daher immer anschließend eine Modelldiagnostik durchführen!

Wie gut ist das Modell?

Zur Erinnerung: das Bestimmtheitsmaß \(R^2\)

Die statistische Kenngröße \(R^2\) ist ein Maß für den Anteil der Variabilität in y, welcher durch das lineare Modell erklärt werden kann:

\[R^2 = \frac{SS_{Regression}}{SS_{Gesamt}} = \frac{\sum(\hat{y}_i-\bar{y})^2}{\sum(y_i-\bar{y})^2}\] \[1 \geq R^2 \geq 0\]

  • Werte nahe 1 weisen auf eine sehr gute Beschreibung der Daten.
  • Werte nahe 0 weisen auf eine sehr hohe Streuung der Daten um die Gerade an
  • Bei \(R^2=1\) lägen alle Datenwerte exakt auf der Geraden.
  • In der Ökologie wird bereits ein \(R^2\) von 0.5 bzw. 50% als sehr gut erachtet.
  • \(R^2=1\) ist auch das Quadrat des Korrelationskoeffizienten \(r\)

Das Bestimmtheitsmaß \(R^2\) im summary() output

summary(mod)
Call:
lm(formula = Petal.Length ~ Sepal.Length, data = iris)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.47747 -0.59072 -0.00668  0.60484  2.49512 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -7.10144    0.50666  -14.02   <2e-16 ***
Sepal.Length  1.85843    0.08586   21.65   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8678 on 148 degrees of freedom
Multiple R-squared:   0.76, Adjusted R-squared:  0.7583 
F-statistic: 468.6 on 1 and 148 DF,  p-value: < 2.2e-16

Modellvalidierung

Annahmen bzw. Voraussetzungen für lineare Regressionsmodelle


  1. Unabhängigkeit (am Wichtigsten!)
  2. Homogenität / homogene Varianz (am Zweitwichtigsten)
  3. Normalität / Normalverteilung
  4. Linearität

Annahmen/Voraussetzungen

Unabhängigkeit

  • Freiheitgrade inkorrekt (Anzahl der unabhängigen Komponenten) → dadurch Beeinflussung der p-Werte (meist stat. Fehler 1. Art erhöht)
  • Gibt
    • Abhängigkeit die man “riechen” kann
    • Abhängigkeit durch nicht passende Modelle
    • Zeitliche oder räumliche Abhängigkeit

Annahmen/Voraussetzungen

Varianzhomogenität

Es wird angenommen, dass

  • die Varianz in Y konstant ist (z.B. die Varianz bleibt unverändert trotz größer/kleiner werdender y-Werte)
  • nur eine Varianz berechnet werden muss und nicht für jeden x-Wert.
  • die Residuen hinsichtlich jeder X-Variable homogen sind.

p

Annahmen/Voraussetzungen

Replikation

Normalität und Varianzhomogenität bei Replikation offensichtlicher

Angenommen Du hast 12 Messwerte für die Variable x und y und willst den (positiven) Einfluss von x quantifizieren:

  • → Überprüfung beider Annahmen anhand einer Exploration der y Variable schwierig
  • → Was aber, wenn es pro X-Wert 50 Beobachtungen gäbe?

Normalität und Varianzhomogenität bei Replikation

Vergleich

Die jeweils 50 Beobachtung bei x = 2,3,4… sollten normal verteilt sein und die Varianzen zwischen den jeweils 50 Beobachtungen homogen:

Normalität und Varianzhomogenität bei Replikation

Vergleich

Die jeweils 50 Beobachtung bei x = 2,3,4… sollten normal verteilt sein und die Varianzen zwischen den jeweils 50 Beobachtungen homogen:

Annahmen/Voraussetzungen

Linearität

  • Wenn eine gerade Linie nicht gut die Beziehung zwischen 2 Variablen beschreiben kann,
    • sollten die Variablen transformiert werden für die Anwendung einer linearen Regression
    • oder es kann eine polynomiale Regression als Spezialfall der multiplen linearen Regression angewendet werden.
  • Die quadratische Funktion z.B. hat drei Parameter und kann verschiedene Typen von Parabeln produzieren: \(y=a+bx+cx^2\)

Wie können die Annahmen überprüft werden?

VOR dem Modellieren

  • Varianzhomogenität:
    • Streudiagramm mit X vs. Y
  • Normalität:
    • Histogramm von Y
    • Shapiro-Wilk-Test mit Y
  • Linearität
    • Streudiagramm mit X vs. Y

NACH dem Modellieren

Führe eine Modellvalidierung durch und erstelle die sog. Diagnostikplots mit den Residuen:

  • Unabhängigkeit:
    • Plotte Residuen vs. jede X-Variable
    • Zeitserien: Plotte die ACF Funktion der Residuen, wende den Durbin-Watson-Test auf die Residuen an
    • Räumliche Daten: Erstelle ein Variogramm mit den Residuen
  • Varianzhomogenität:
    • Streudiagramm mit Residuen vs. vorhergesagte Y-Werte
  • Normalität:
    • Histogramm oder QQ-Plot der Residuen

Überprüfe Ausreißer/einflussreiche Beobachtungen

Leverage / Cook’s Distanz ⇒ weiterer Schritt der Modellvalidierung

  • ‘Leverage’ h: ein Werkzeug, dass Extremwerte der erklärenden Variable identifiziert, welche potentiell auch das Regressionsergebnis beeinflussen könnten.
  • Cook’s Distanz D: das Maß für zu einflussreiche Datenpunkte. Identifiziert einzelne Beobachtungen, die alle Regressionsparameter beeinflussen. Es berechnet eine Teststatistik D, die die Veränderung aller Regressionsparameter misst, wenn eine Beobachtung ausgeschlossen wird.
    • D > 0.5 kritisch
    • D > 1.0 sehr kritisch

→ Es ist einfacher den Ausschluss beeinflussender Datenunkte zu rechtfertigen, wenn eine große Cook’s Distanz und eine große Leverage für den Datenpunkt vorliegt.

Überprüfe Ausreißer/einflussreiche Beobachtungen

Leverage / Cook’s Distanz ⇒ weiterer Schritt der Modellvalidierung


p

Standardgrafiken zur Modellvalidierung

par(mfrow = c(2,2))
plot(mod)

Gewöhnliche vs standardisierte Residuen

  • Gewöhnliche Residuen (\(e_i\)): beobachtete Werte - vorhergesagte Werte (‘observed – fitted’)
  • Standardisierte Residuen (\(e_{stand.}\)): Residuum dividiert durch seine Standardabweichung: \[e_{stand.} = \frac{e_{i}}{\sqrt{MS_{Residual}*(1-h_{i})}}\]
    • wobei \(h_{i}\) = leverage der einzelnen Beobachtungen i; \(MS_{Residual}\) = Mittlere Summenquadrate der Residuen
  • Standardisierte Residuen werden als normalverteilt angenommen mit einem Mittelwert von 0 und einer Varianz von 1; N(0,1).
  • Folglich: große Residuenwerte (>2) bedeuten eine schlechte Anpassung des Regressionsmodels.

Berechnen der Residuen in R

Du kannst beide Arten von Residuen berechnen mit

  • residuals(model) (oder resid()) vom ‘stats’ Paket → ergibt einen Vektor mit den gewöhnlichen Residuen

  • rstandard(model) vom ‘stats’ Paket → ergibt einen Vektor mit den standardisierten Residuen

  • add_residuals(data, model, var = "resid") vom ‘modelr’ Paket (in ‘tidyverse’) → fügt die Variable ‘resid’ mit den gewöhnlichen Residuen zu deinem ‘data frame’ hinzu; hilfreich beim ‘piping’ von Funktionsaufrufen!

Berechnen der vorhergesagten Werte in R

Du kannst Dir die vorhergesagten Werte für dein Modell über 3 Funktionen ausgeben lassen:

  • fit(model) → gibt einen Vektor mit den angepassten Werte zurück, die vom Modell für die Werte der erklärenden Variablen vorhergesagt werden
  • predict(model) von dem internen ‘stats’ Paket → verwendet Informationen aus dem angepassten Modell, um eine Glättungsfunktion für die Darstellung einer Linie im Streudiagramm zu erstellen. Der Funktion kann ein data frame mit Sequenzen von Werten aller X Vaiablen, die im Modell sind, übergeben werden (Argument newdata).
  • add_predictions(data, model, var = "pred") von dem ‘modelr’ Paket (in ‘tidyverse’) → fügt die Variable ‘pred’ mit den vorhergesagten Werten zu deinem Datensatz hinzu; hilfreich bei ‘piping’ Operationen!

Weitere Diagramme zur Modellvalidierung

Modellanpassung

Vorhergesagte vs. beobachtete Werte

iris %>% modelr::add_predictions(mod) %>%
  ggplot(aes(x = pred, y = Petal.Length)) + 
  geom_point() + geom_abline()

Weitere Diagramme zur Modellvalidierung

Residuen vs. jedes X

iris2 <- iris %>% modelr::add_residuals(mod)
p <- ggplot() + geom_hline(yintercept = 0, color = "red")
p1 <- p + geom_point(aes(x = Sepal.Length, y = resid), iris2)
p2 <- p + geom_point(aes(x = Sepal.Width, y = resid), iris2)
p3 <- p + geom_point(aes(x = Petal.Width, y = resid), iris2)
gridExtra::grid.arrange(p1,p2,p3, nrow = 1)

→ Die Residuen sollten zufälligt verstreut sein. Ein Muster deutet auf eine Fehlspezifikation des Modells hin, wie hier zum Beispiel.

Visualisierung der linearen Regression

mit ggplot2

Verwende geom_smooth und setze das Methodenargument als Name der Modellierungsfunktion (hier: ‘lm’). Falls Du den Standardfehler nicht anzeigen willst, setze zusätzlich se = FALSE:

iris %>% ggplot(aes(x = Sepal.Length, y = Petal.Length)) +
  geom_point() + geom_smooth(method = "lm", se = FALSE)

Visualisierung der linearen Regression

manuell

iris %>% modelr::add_predictions(mod) %>%
  ggplot(aes(x = Sepal.Length)) + geom_point(aes(y = Petal.Length)) + 
  geom_line(aes(y = pred)) 

Your turn …

Shiny App

Quiz: Verteilung & Modellspezifikation

Erkunde die Shiny App auf der folgenden Folie und beantworte diese zwei Fragen:

Shiny App (scrolle weiter runter)

Fragen?